Cloud, Shadow and Snow Masking using Multispectral Signature Analysis

There are many reasons for identifying and masking clouds. This technique is used to assess the quality of archived images, for example Landsat GLOVIS. Sometimes is necessary to flag cloud and cloud shadow areas for land use analysis so that cloudy or cloud shadow areas will not be included in subsequent image processing. Another interesting idea is to synthesize pixels in the cloudy areas. This can sometimes be done as the cloudy areas contain less but not zero information. It is also often necessary to discriminate between clouds and snow cover in making such assessments.

The standard PANCROMA™ utilities for masking clouds use either three or five of the Landsat multispectral band files. These files are used in a non-linear combinatorial manner to predict which image pixels correspond to cloud cover and which do not. These utilities, particularly the five-file utility based on the Luo method can be highly effective.

However, these techniques are only effective in targeting clouds. An alternative approach to creating masks is to use the PANCROMA™ multispectral analysis tools. The advantage is that virtually any type of land feature can be targeted (for example cloud shadows and snow) using a straightforward technique. These tools can also be highly discriminating. This paper describes how this can be done. In the case described in this article, cloud, cloud shadow and snowy areas will be masked out of an image so that pixels from another image can be substituted in order to improve its appearance. In the process, it will illustrate how a wide variety of surface features can be identified using PANCROMA™ spectral analysis tools.

Let's say that you have a Landsat image that interests you but that image is obstructed with cloud and snow cover. One alternative is to abandon the image and look for one collected earlier or later than your target image. Another alternative is to use the image that you favor, but substitute pixels from another image to improve the appearance of your target. The idea is to substitute as few pixels as is necessary to improve the appearance of the target while maintaining as much of the original character of your target image as possible.

The first step in this procedure is to download all of the multispectral bands for both your target image, and for the image from which you will 'borrow' pixels. I will call this second image the Reference image. You will need all the bands 1,2,3,4,5,6 and 7 for each image. Note that for many Landsat multispectral bundles, the thermal infrared (TIR) band6 image will be half resolution (approximately 60m/pixel) and half the size of the multispectral bands. Some band6 data has been processed by expanding the image and interpolating the gaps. If your data is archived as a 60m image, you can use the PANCROMA™ TIR interpolation utility to expand and interpolate it to match the size and scale of the multispectral files. See the PANCROMA™ Instruction Manual for directions on how to do this. My TIR bands had already been processed by the USGS to this resolution.

For my example I chose two mountain desert images from northern Chile. My target image is obstructed with clouds, cloud shadow and snow. I have another image that has zero cloud cover and much less snow. I will substitute pixels from the clean image into my target. These two images are shown below. I chose my target image to be considerably obscured by cloud and snow in order to better illustrate the technique. As you can see, the images are subsets of the Landsat scene. I used the PANCROMA™ utilities; subsetting to consistent corner coordinates three band files at a time. Instructions for doing this are also given in the PANCROMA™ Instruction Manual.

Target Image obscured with cloud, cloud shadow and snow

Reference Image - we will "borrow" pixels from this image

We will first generate spectral signatures using the PANCROMA™ Point Spectrum Generator. We will do this three times, once for the clouds, once for the cloud shadows, and once for the snow. We will then use the generated image spectra and the Euclidean Distance Analyzer to identify our targets for removal, with the masking option switched on.

Begin by opening your seven Landsat Target bands in the order: band1, band2, band3, band4, band5, band6 and band7 by selecting 'File' | 'Open' from the main menu. Now select 'Spectral Analysis' | 'Landsat Point Spectrum Generator' | 'Six/Seven 8-Bit DN Bands'. The Band Display Selection Form will become visible as shown below. This form allows you to specify which band or band combination you want to use to identify your target pixels. We will be going after the clouds first, so band1 is a good choice, although others will work as well. Click 'OK' when you have selected your band.

Band Selection Form - choose one of the grayscale bands or a three-band composite for pixel selection. I selected band1 for this example because Rayliegh scattering makes clouds and atmospheric water vapor show up better than in other bands.

Now the selected band image will appear, along with the Point Spectrum Form. The main purpose of the Point Spectrum Form is to keep a running tally of your selected points. It also controls the zoom feature in case you need to get in closer into your image to accurately identify your target pixels. Using your cursor, successively click on the cloudy areas of the image until you have collected a representative sample (You can choose as many as 50, but you can often get by with far fewer.) I chose 18. When you are finished, select 'OK'.

Point Spectrum Form - you will click on the band image to identify target cloud pixels. The TOA reflectance of each pixel will be computed and used to generate a target spectrum.

Since we will be computing a Top of Atmosphere (TOA) Reflectance spectrum, the TOA Reflectance Data Form will become visible. Enter the Solar Elevation Angle, Julian Day and TIR Gain Level corresponding to your image from the Landsat metadata file. I left the 'Landsat ETM+' radio button checked since this is my image data type. Check the 'Retain Settings' box so that the entered values will be saved for subsequent runs. See the Instruction Manual for details on how to find the metadata numbers.

TOA Form - in order to compute TOA reflectances from image digital numbers (DNs) metadata parameters Solar Elevation Angle and Julian Day of acquisition are required. Other data needed to compute TOA reflectances like the Gain and Bias values are always the same for ETM+ bands so can be stored in the program.

Now select 'OK' and the reflectance spectrum for your targeted pixels will be generated as shown below. This sample will be used to identify all pixels of that type (cloud pixels in this case) in the Target image.

Target Spectrum - a plot of reflectances versus band number for the target pixels that you selected. Average reflectances are shown in red. This spectrum will be used to identify all of the cloud pixels in the data bundle that meet our selection criteria.

The next step is to generate the first mask. Select 'Close Graphics Window and Reset' from the Main Menu. When you do so, the spectrum will close. However the average TOA reflectances will be stored in the Euclidean Distance Analyzer for the next step. Re-open your seven Target band files as before by selecting 'File' | 'Open' from the main menu. Select 'Spectral Analysis' | 'Landsat Spectral Analyzer' | 'Six/Seven 8-Bit DN Files' | 'Euclidean Distance'. (Note, the Manual Method can also be used for this step. We will illustrate mask generation using the Euclidean Distance Analyzer in this example.) The Spectral Criteria form will become visible as shown below.

Spectral Criteria Form - the target vector (based on the target spectrum above) has been automatically entered into the form. Just check the 'Create Mask' check box and set the Level Slider.

Note that there are seven track bars corresponding to the seven multispectral bands we have opened. The average reflectances are already input into these track bars so no adjustment is required (unless you feel some is necessary). Check the 'Retain Settings Between Runs' check box. This will hold our Spectral Criteria settings in case we need multiple runs to make our mask. Also check the 'Create Mask' check box in the Mask Group. When you do so, the Level Selection track bar in the Mask Group will become enabled.

This track bar lets you specify how far away the Euclidean Distance vector (which will be computed for each pixel in the multispectral data set) can be from the target vector (as defined by the band Track Bar settings) and still be included in the mask. The higher the number from 0 to 8, the farther the allowable distance and the more pixels will be included in the mask. Low numbers mean high discrimination but possibly too narrow a specification with some valid pixels being excluded. High number mean lower discrimination with the possibility of including unwanted pixels in your mask. I will use Level=3 for my first run. Select 'OK'.

When you do so, the TOA Reflectance Data form will again become visible. The values you used for the Point Spectrum Generator and click 'OK' and the mask will be generated. The color contour plot and the mask show the Euclidean distances between each pixel reflectance and the spectral signature target reflectance. The mask image has zeroed out image pixel values at locations corresponding to the cutoff you specified with the Level Selection track bar.

Cloud Mask - the best match pixels (closest Euclidean Distance) are red. We have included up to the color yellow in the mask by setting the Level Track Bar onto position '3'.

The analysis has been very successful at discriminating the clouds from everything else in the image. In particular, note how the bright cloudy areas are identified, while the bright snowy areas are not included. To the human eye, these areas are indistinguishable. To the multispectral "eye" of the Landsat sensors, they are quite different. However the analysis has not picked up the cloud "fringe" areas very well where the clouds are less dense and where there is a transition to clear atmosphere. We can consider alternative approaches to engaging the transition pixels without selecting the snow pixels in the process. I will discuss this in a future article. For this exercise I simply increased the Level Track Bar a couple of notches until I got them all, causing some overlap into the snowy areas. Since I needed to mask these anyway this was not a problem for this image.

The next step is to save the mask. Select 'File' | 'Save Grayscale Image' | 'GeoTiff' from the menu and give the file a name such as 'cloudMask.tif'. Select 'Close Graphics Window and Reset'. Now repeat the steps outlined above for the cloud shadow areas. This is done exactly as described except rather than clicking on the cloudy areas, you will use the Point Spectrum Generator to click on the cloud shadow areas. Then create a cloud shadow mask and save it as 'shadowMask.tif' or the like. Note that in addition to the cloud shadows, I have picked up some shadows cast by the mountain peaks. Since a shadow is a shadow, this is to be expected and I will accept these unwanted but harmless additions. Note that I was highly selective with my Level Track Bar setting for the cloud shadows, moving it to level '0' so I would only pick up the shadows.

Cloud Shadow Mask - a little hard to see but the shadows are in bright red. You have to look closely to see the zero mask values in the grayscale image but they are there.

Now repeat the steps in the preceding paragraph for the snowy areas, clicking on these in the Point Spectrum generator to generate a snow spectrum and then using the Euclidean Distance Analyzer to create a snow mask. Save this grayscale mask as 'snowMask.tif'. My snow mask is shown below.

Snow Mask - the snowy areas have been identified without being confounded with the clouds. Need a little more image processing programming to get rid of the fringes around the snow fields.

The next step is to create a composite mask, combining the masked areas in the cloud mask, the shadow mask and the snow mask into a single mask, which we will use to blank out the unwanted regions in our target bands. To do this, open the cloud mask and the shadow mask by selecting 'File' | 'Open' and then 'Preprocess' | 'Band Math' | 'Boolean AND Two Images'. A new image will be created that will add the zero value pixels of the shadow mask to those of the cloud mask. Save this intermediate mask as 'cloudShadowMask.tif'.

Select 'Close Graphics Window and Reset' and then open the cloudShadowMask.tif file and snowMask.tif by selecting 'File' | 'Open' and then 'Preprocess' | 'Band Math' | 'Boolean AND Two Images' again. At the end of this operation, you will have your composite mask, which you can save as 'cloudShadowSnowMask.tif. My composite mask is shown below.

Composite Mask - created by Boolean 'AND'ing the three individual masks together.

Now we will mask our target images. Select 'Close Graphics Window and Reset' and open your first target band file (i.e. the one with the clouds and snow and the 'cloudShadowSnowMask.tif by selecting 'File' | 'Open'. Boolean AND the mask and the band file image by selecting 'Preprocess' | 'Band Math' | 'Boolean AND Two Images'. Save the file as band1MaskedImage.tif. . Select 'Close Graphics Window and Reset' and repeat for your band2 and band3 target images, saving each masked image in turn. You should now have three images, each with the cloud, shadow and snowy areas masked out as shown in the image below.

The final step is to fill the gaps. Open the band1MaskedImage.tif file and your band1 Reference (clean) image using 'File' | 'Open'. The fill the gaps by selecting 'Gap Fill' | 'Gap Fill Transfer Method'. PANCROMA™ will transfer pixels from your clean image into the zero pixel-masked gaps in your target image. Save the band file as 'filled1.tif'. Select 'Close Graphics Window and Reset' and repeat for target band2 and band3. The image shown below is the result of my filling my band1 target image. A histogram spectrum is also shown as a global histogram match is normally conducted in order to get the image tones between the Target and Reference fill pixels closer.

Filled Image - used the simple PANCROMA™ Transfer method: fill in the gaps from the Reference image and conduct a global histogram match so the color tones are close. We could have used a more sophisticated algorithm but this was an easy match and the Transfer method is fast.

I created an RGB color composite of my processed image, which is shown below. Compare it with the first image in the article. While not perfect, it gives you some idea of what you might accomplish with images that are not so extensively obscured as this one. It also illustrates the power of PANCROM™ multispectral analysis tools for discriminating various types of ground targets other than clouds and snow. PANCROMA™ has similar utilities for masking Landsat Reflectance, SPOT^®, ASTER, EO-1 ALI, and WorldView-2^® data.

Processed Image - the fruits of our labors. Needs additional image processing to clean it up. To be continued...